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NOMENCLATURE 


* 

a 

r 

reference  speed  of  sound  (cm/sec) 

c 

Pi 

"k  k k*? 

c T:  /a  I specific  heat  at  constant  pressure 

p r r 

°12 

D, „ /a  L ; coefficient  of  diffusion  (used  in  Pick's  Law) 

12  r 

E 

* , *2  u^  / 

E /a^  ; total  energy  = n + 2 P'P 

E 

Activation  energy  (kcal/mole) 

h 

h /a  ; enthalpy,  e Y h 

^ i=l 

h. 

1 

I c dT  + h.  ; enthalpy  of  species  i (thermal  + chemical) 

JjO  Pi 

h.° 

1 

heat  of  formation  of  species  i at  T 

* 

AH 

heat  of  reaction  (kcal/mole),  h^^  - for  A 2B 

k 

8 

k T /p  a L ; thermal  conductivity 

g r '"^r  r 

L 

reference  length  (cm) 

M. 

1 

molecular  weight  (gm/gm-mole) 

N 

total  number  of  species 

P 

k k k^ 

P /Pj.  \ ; pressure 

%fall 

q*  a ; external  heat  flux,  q n “ cal/cm  -sec 

^ wall  ^r  r wall 

S^all 

o 

constant  portion  of  external  heat  flux 

Re 

k k k k 

p^  L /y^  ; reference  value  of  Reynolds  number 

R 

o 

k k k2 

R T /a  ; universal  gas  constant 

0 r r 

t 

k k k 

t a /L  : time 
r 

T 

k k 

T /T^  ; temperature 

* 

T 

r 

reference  value  of  temperature  (°K) 
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NOMENCLATURE  (Continued) 


u 


Y. 

1 


ic  * 

u /a^  ; velocity 

p^/p;  mass  fraction  of  species  i 


R 


Pi 


* ^ 

z /L  ; axial  coordinate 
pre-exponential  factor  (sec 

"k  k 

y /y^  ; coefficient  of  viscosity 
(4/3)vi/Re 

reference  value  of  coefficient  of  viscosity  (gm/cm-sec) 

k k 

p /p^  ; density 
density  of  species  i 

3 

reference  value  of  density  (gm/cm  ) 


a constant,  0 0 ^ 1 

k 

~i  * 

0)  L /p  a ; chemical  production  rate  of  species  i 
i 


Subscripts 


i pertaining  to  species  i 

j location  of  grid  point 

r reference  value 

t derivative  with  respect  to  time 

z derivative  with  respect  to  axial  coordinate 


Superscripts 

^ dimensional  quantity 

n current  time  level 

n+1  new  time  level 
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I.  INTRODUCTION  AND  MOTIVATION 


Any  attempt  to  construct  a mathematical  model  of  solid  propellant 
ignition  and/or  unsteady  combustion  is  faced  with  the  difficulty  of  des- 
cribing the  gas-phase  reaction  zone.  Unless  it  can  be  assumed  that 
solid-phase  reactions  dominate  completely,  a prediction  of  solid  pro- 
pellant decomposition  rate  is  heavily  dependent  upon  events  occurring 
in  the  gas  phase.  Competition  between  convection,  heat  conduction, 
species  diffusion,  and  chemical  reaction  rates  will  establish  the  in- 
stantaneous energy  flux  to  the  surface  of  the  reactive  solid.  The 
ability  to  describe  this  interaction  is  a prerequisite  for  mathematical 
models  designed  for  transient  calculations. 

The  study  of  transient  solid  propellant  combustion  in  Ref.  1 ex- 
amined several  combustion  models  which  assume  the  flame  region  to  be 
quasi-steady.  This  affords  a tremendous  simplification  in  the  analysis, 
if  it  also  can  be  assumed  that  the  spatial  distributions  of  heat 
release  and  chemical  production  rate  are  known  and  remain  unchanged. 

The  latter  assumption  effectively  uncouples  the  energy  equation  from 
the  species  continuity  equations,  and  thus  any  fixed  distribution  is 
hard  to  defend.  Furthermore,  the  basic  quasi-steady  assumption  is 
difficult  to  justify  for  problems  related  to  a gun  combustion  chamber 
where  the  entire  ignition  and  unsteady  combustion  process  is  completed 
within  milliseconds. 

2 3 4 5 ... 

Four  recent  models  * * * of  the  solid  propellant  ignition  process 

have  analyzed  the  gas-phase  reaction  zone  in  much  greater  detail.  These 
studies  have  led  to  considerable  insight  into  the  competing  processes 
which  control  ignition  of  a reactive  solid.  However,  all  have  dealt 
with  a constant-pressure  "isolated"  ignition  problem.  The  isolation  is 
a consequence  of  applying  vanishing  gradient  boundary  conditions  to  the 


Kooker^  D.  , and  Nelson^  C.  IV.,  ^'Numerioal  Solution  of  Three  Solid 
Fvo'pellant  Combustion  Models  During  a Gun  Pressure  Transient^^j  BRL 
Report  No.  1955  (AD  #A0Z5250)^  USA  Ballistic  Rsoh  Lab^  APG,  MD^ 

Jan  77  (see  also  12th  JANNAF  Combustion  Meeting^  Newport ^ RI^  CPIA 
Publication  273j  pp  175-197 ^ Dec  75). 

^ Hermanoej  C.  W.  j and  Kumar,  R.  K. , "Gas  Phase  Ignition  Theory  for 
Homogeneous  ProipetZants  Under  Shook  Tube  Conditions",  AIAA_J_. , Vot. 

8,  No.  9,  pp  1551-1558  (Sep  70). 

^ Kashiwagi,  T.,  "A  Radiative  Ignition  Model  of  a Solid  Fuel",  Comb. 
Sci.  Tech.,  Vol.  8,  pp  225-236  (1974). 

^ Bradley,  H.  H.,  Jr.,  "A  Unified  Theory  of  Solid  Propellant  Ignition, 
Part  1,  Development  of  Mathematical  Model",  NWC  TP5618  (Aug  74). 

^ Kindelan,  M.,  and  Williams,  F,  A.,  "Radiant  Ignition  of  a Combustible 
Solid  with  Gas-Phase  Exothermiaity",  Acta  Astronautioa,  Vol.  2, 

No.  11/12  Nov /Deo  75). 
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gas-phase  dependent  variables  at  an  infinite  distance  from  the  pro- 
pellant surface;  thus,  changes  in  the  outside  environment  cannot  in- 
fluence the  ignition  process.  Additionally,  each  analysis  is  constrained 
by  various  combinations  of  assumptions  such  as  a unity  Lewis  number, 
negligible  kinetic  energy,  multiple  species  with  identical  molecular 
weights  and  specific  heats,  etc. 

The  current  study  is  motivated  by  a desire  to  model  the  ignition 
and  unsteady  combustion  process  of  an  end-burning  solid  propellant  grain  ' 

confined  in  a one-dimensional  chamber,  as  shown  in  Fig.  la.  The  model 
must  be  sufficiently  general  to  allow  variable  transport  properties, 
various  chemical  decomposition  schemes,  etc.,  without  the  restriction  v 

of  constant  pressure.  The  objective  is  two-fold:  (a)  predict  the  in- 
stantaneous regression  rate  of  the  reactive  solid  as  a function  of  the 
self-generated  pressure,  and  (b)  look  for  potential  instabilities  in 
both  regimes — ignition  and  unsteady  combustion. 


(a) 


^wall 


(t) 


REACTIVE  GAS 


(b) 


Figure  1.  Problem  Configurations  - (a)  Solid  Propellant  in  Closed 
Chamber;  (b)  Reactive  Gas  in  Closed  Chamber 

The  present  investigation  is  the  first  step  in  the  above  model 
and  has  the  limited  objective  of  describing  the  ignition  sequence  of  a 
reactive  gas  confined  in  a rigid,  one-dimensional  chamber.  As  shown  in 
Fig.  lb,  one  boundary  is  assumed  to  be  insulated  and  the  other  is  sub- 
jected to  an  externally-prescribed,  time-dependent  heat  flux.  The  im- 
posed heat  flux,  as  used  here,  serves  as  the  ignition  stimulus  and  is 
not  intended  to  represent  the  response  of  the  solid  in  the  full  pro- 
pellant ignition  problem  of  Fig.  la.  The  remainder  of  this  report  is  t 

concerned  with  the  solution  of  the  gas-phase  ignition  problem  shown  in 
Fig.  lb. 


II.  ANALYSIS 

The  gas-phase  ignition  problem  described  above  requires  the  numeri- 
cal solution  to  the  one-dimensional,  time-dependent,  compressible 
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Navier-Stokes  Equations  for  a multi-species  reacting  flow  field.  To 
correctly  model  the  influence  of  the  finite- size  chamber,  the  system 
must  retain  the  full  momentum  equation,  as  opposed  to  the  trivial 
constant-pressure  form.  At  no  time  will  the  analysis  employ  an  externally 
specified  value  of  "T.  thus,  no  artificial  boundary  is  drawn 

between  ignition  and  unsteady  combustion.  Ignition  will  be  determined 
by  a go/no-go  response  when  the  external  heat  flux  is  terminated. 

Assumptions 

The  reactive  gas  confined  in  the  chamber  is  assumed  to  be  a mixture 
of  perfect  gases,  each  of  which  may  have  a different  molecular  weight, 
specific  heat  capacity,  and  heat  of  formation.  It  is  assumed  that  the 
mixture  viscosity,  thermal  conductivity,  and  species  diffusion  co- 
efficients can  be  described  at  the  molecular  level,  i.e.,  all  transport 
is  laminar.  Admittedly,  turbulent  transport  is  the  dominant  mode  in 
many  combustion  situations.  However,  the  gas-phase  ignition  problem 
posed  here  should  not  be  controlled  by  turbulence.  In  the  event  of  a 
violent  unsteady  flame  propagation  in  the  chamber,  this  assumption 
would  have  to  be  re-examined.  The  numerical  solution  procedure  assumes 
variable  laminar  transport  properties,  with  no  constraints  on  the  magni- 
tude or  variability  of  the  mixture  Lewis  number  or  Prandtl  number.  The 
influence  of  radiation,  gravity,  bulk  viscosity,  the  Dufour  effect,  the 
Soret  effect,  and  diffusion  due  to  pressure  gradients  will  be  neglected. 

The  current  results  are  based  on  the  additional  assumption  that  the 
chamber  initially  contains  a reactive  gas  composed  of  the  single  species 
A,  capable  of  undergoing  the  one-step,  irreversible  reaction 

A 2B. 

Pick's  Law  is  then  used  to  describe  the  diffusion  mass  flux  in  this  bi- 
nary system.  It  should  be  noted  that  extending  Pick's  Law  to  a generalized 
multi-component  mixture,  and  then  assigning  unequal  diffusion  coefficients, 
will  lead  to  a violation  of  the  basic  identity  that  the  sinii  of  the  dif- 
fusion mass  fluxes  must  be  zero.  Although  this  has  been  done^>^,  the 


Rivard^  W.  C. , Farmer j 0.  A.,  and  Butter ^ T.  D.^  "RICE:  A Computer 
Program  for  Multicomponent  Chemically  Reactive  Flows  at  All  Speeds''^ 
LA-58123  Los  Alamos  Scientific  Laboratory,  Los  Alamos,  EM,  Mar  75. 
y 

Kothari,  A.  P.,  Anderson,  J.  D.,  Jr.,  and  Jones,  E.,  "Navier-Stokes 
Solutions  for  Chemical  Laser  Flows",  AIAA  J.,  Vol.  15,  No.  1, 
pp  92-100  (Jan  77). 
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correct  procedure  requires  the  solution  of  the  coupled  equations  relating 
all  diffusion  mass  fluxes  and  mass  fraction  gradients,  as  described  by 
Edelman^ . 

I 

j 

Governing  Equations  ' 

Based  on  the  foregoing  assumptions,  the  conservation  equations  in 
an  inertial  frame  of  reference  can  be  written  as: 


Global  Continuity 

P^.  + (pu)^,  = 0 

(1) 

Momentum 

(2) 

Energy 

(pE)^  + [ojE  + pu]^  = + p^uu^ 

^ ) 

+ >}z 

z 

(3) 

Species  Continuity 

pY.  + (i)Y.  = (pD,_Y.  ) + 0)^ 

1 1 12  1 z 

t z z , - 

(4) 

i=i,. 

N 

. , ,N- 

State 

P = pRT  , Re  R .E, (Y./M.) 

O 1=1  1 1 

N r T '' 

h = .E.Y.h.  , h.  = / ^ c dT  + h.° 

1=1  11  " 1 / ^ p.  1 

2 i 

(5) 

where  to  = pu,  E = h + ^ - p/p,  = P^/p 


Equations  (1)  - (3)  are  solved  numerically  in  their  strict  conservation 
form.  Attempts  to  solve  the  species  continuity  equation  written  in 
terms  of  partial  densities  led  to  numerical  instabilities  when  forming 
a mass  fraction  with  a value  very  close  to  unity.  The  origin  of  this 
instability  appears  to  be  machine  inaccuracy  in  division  of  two  nearly 
identical  numbers,  and  is  not  attributable  to  the  stiffness  of  the 
equation.  Use  of  the  form  shown  in  Eq.  (4)  eliminated  the  difficulty. 


The  common  practice  of  transforming  the  axial  coordinate  (using 
an  integral  of  density)  to  eliminate  direct  dependence  on  the  global 
continuity  equation  is  unworkable  here  because  the  computational  domain 
is  enclosed  by  two  boundaries  at  finite  locations. 


Edelman^  R.  B.^  Fovtuney  0.  and  Weilerstein^  ”An  Analytical 
and  Eccpevimental  Investigation  of  Gravity  Effects  Upon  Lconinar  Gas 
Jet-Diffusion  Flames”^  14th  Symposiim  (International)  on  Combustion^ 
Aug  72j  Penn  State  University ^ pp  399-412, 
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Boundary  Conditions 


Boundary  values  of  the  dependent  variables  must  be  specified  at  both 
the  heated  boundary,  z=0,  and  the  insulated  boundary,  z=l  (see  Fig.  2). 


■ fci 


z 

(a)  u H 0 u E 0 

(b)  = 0 i=l,,..N-l 

<<=>  "I  ■ 

Figure  2.  Physical  Boundary  Conditions 

The  assumption  of  impenetrable  walls  requires  the  mean  flow  velocity  to 
vanish  at  z=0  and  z=l  [Eq.  C6a)].  The  individual  species  velocities 
must  also  vanish  at  the  solid  boundary,  which  means  that  all  diffusion 
velocities  are  zero.  Since  by  Fick’s  Law  the  diffusion  velocity  is 
linearly  dependent  upon  the  spatial  gradient  of  the  mass  fraction,  each 
mass  fraction  gradient  must  vanish  at  z=0  and  2=1  [Eq.  (6b)].  Combining 
Eqs.  (6a)  and  (6b)  with  the  energy  balance  at  the  boundaries  results  in 
a zero  temperature  gradient  at  z=l  and  a time-dependent  temperature 
gradient  at  z=0  governed  by  ^vall^^^  (6c)].  At  each  boundary, 

there  are  N+3  dependent  variables,  but  only  N+2  conditions  given  by 
Eq.  (6)  along  with  the  equation  of  state.  However,  no  physical  boun- 
dary condition  is  missing.  The  situation  is  directly  analogous  to  the 
simpler,  inviscid  acoustic  problem  at  a boundary;  a value  of  the  velo- 
city perturbation  (or  an  equivalent  admittance  condition)  can  be  speci- 
fied but  the  value  of  the  pressure  perturbation  must  follow  from  the 
solution.  To  specify  both  values  would  over-determine  the  system  of 
equations . 

The  numerical  solution  procedure  requires  values  of  all  dependent 
variables  at  both  boundaries  for  every  time  step.  In  a complex  system 
of  equations,  the  manner  in  which  the  "extraneous”  boundary  values  are 
determined  is  usually  the  difference  between  success  and  failure.  An 
accurate  method  to  solve  the  inviscid  acoustic  problem  is  based  on  the 
fact  that  the  system  of  equations  is  hyperbolic.  Because  of  the  finite 
propagation  velocity,  one  characteristic  line  in  the  t-z  plane  will  in- 
tersect the  boundary  at  the  new  time  level.  The  compatibility  condition 
along  this  line  provides  the  missing  information  which  relates  the  un- 
known boundary  values  to  the  known  solution  at  the  previous  time  level. 


Y.  = 0 (6) 

z 

T = 0 
z 
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In  the  present  gas-phase  ignition  problem,  the  governing  equations  are 
parabolic.  All  characteristic  directions  have  essentially  collapsed 
to  the  single  line  t=constant,  which  implies  an  infinite  speed  of 
propagation.  The  potential  benefit  from  a method-of-characteristics 
solution  at  the  boundary  is  now  lost. 


The  present  investigation  attempted  several  different  methods  to 
determine  the  boundary  value  of  pressure  consistent  with  the  interior 
numerical  solution.  The  only  successful  one  is  based  on  the  idea  of 
mechanical  equilibrium  at  the  boundary.  Since  the  mean  velocity  is  con- 
strained to  be  zero  at  all  times,  the  momentum  equation  reduces  to  the 
statement  that  the  gradient  of  the  stress  tensor  must  vanish.  Thus,  at 
the  boundary,  the  gradient  of  thermodynamic  pressure  must  balance  the 
gradient  of  the  viscous  normal  stress,  viz., 


(7) 


Equation  7 is  enforced  at  a boundary  by  fitting  a cubic  profile  (in 
velocity  and  pressure)  through  the  adjacent  nearest  neighbors,  and 
appropriately  differentiating.  This  provides  the  value  of  pressure, 
and  then  the  value  of  density  follows  directly  from  the  equation  of 
state  [Eq.  (5)]. 

The  boundary  conditions  are  completed  by  specifying  the  time— de- 
pendent heat  flux  imposed  at  z=0.  All  computations  discussed  in  the 
present  study  are  based  on  qwall(t)  shown  in  Fig.  3,  where  the  duration 
of  heating  is  altered  by  varying  the  cut-off  time,  t^q.  The  shut-down 
portion  of  the  heating  profile  (after  t^q)  is  necessary  to  keep  a dis- 
continuous gradient  from  destroying  the  numerical  solution;  its  shape 
remained  unchanged  throughout  the  study.  Unless  otherwise  stated,  the 
constant  portion  of  the  heating  curve  was  held  at  the  value  10  cal/cm^-sec. 
Note  that  for  t > ^^q  5,  the  heated  wall  becomes  adiabatic. 


Figure  3,  Time-Dependent  Heat  Flux  Imposed  at  z=0 
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Method  of  Solution 


The  numerical  integration  technique  used  to  solve  Eqs.  (1)  - (4) 
is  constructed  on  the  basis  of  a control  volume  of  length  Az  centered 
about  each  grid  point  (see  Fig.  4). 


j 1 

I I 

9 ^ 


J-i  I j I 

\m»mnnnmmm7TT^ 


AZ' 


J+1 


Figure  4.  Control  Volume  Nomenclature 

Each  term  in  the  equation,  other  than  the  time  derivative,  is  weighted 
between  its  value  at  the  old  time  level  ”n"  and  the  new  time  level  ”n+l”. 
Using  global  continuity  [Eq.  (1)]  as  an  example. 


!i 


n+1 


n 


At 


+ 


Az 


(pu 


3+4  - ‘’“j-4> 


n+1  . 1-0 


+ 


Az 


'"“3+4 


pu,  - 0 

J 2 


(8) 


where  0 is  a constant,  0 < 0 < 1,  and  the  subscripts  indicate  grid  point 
location.  The  accuracy  of  the  integration  is  O(At^)  if  0=2.  However, 
with  the  competing  effects  of  convection  and  diffusion,  stability  can  be 
assured  only  when  0 > i.  All  solutions  discussed  in  this  report  are 
based  on  0=0.55.  Values  of  the  dependent  variables  are  stored  only  at 
the  integer  grid  point  location;  values  required  at  the  midpoints  are 
found  from 


pu._^  = §(pu  + pu  ),  etc. 

J 


(9) 


Similarly,  flux  terms  which  depend  on  a spatial  gradient  are  represented 
as 


k T 
g z 


= 4(k  +k  )(• 
j+i  ®j+l 


- T. 

2±1 1) 

Az  ^ 


(10) 


using  the  internal  heat  flux  as  an  example.  This  procedure  insures 
strict  conservation,  i.e.,  the  flux  value  entering  control  volume  j+1 
across  the  interface  located  at  j+§  will  be  exactly  the  value  which 
left  control  volume  j across  the  same  interface. 
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Equations  (2)  - (4)  are  each  written  in  a form  analogous  to  Eq.  (8), 
such  that  the  terms  representing  time  storage,  convection,  diffusion,  and 
production  are  coupled  at  the  new  time  level  ”nH-l".  The  result,  after 
applying  the  boundary  conditions,  is  a tri-diagonal  matrix  which  can  be 
inverted  with  the  standard  Thomas  Algorithm.  The  solution  then  proceeds 
by  successive  substitution,  not  quasi-linearization.  A series  of  tri- 
diagonal matrix  inversions  is  performed  in  the  following  sequence,  where 
each  step  incorporates  the  latest  dependent  variable  estimates  from  the 
previous  steps: 

1.  Solve  the  species  continuity  equations  for  the  Y^*s, 

2.  Solve  the  energy  equation  for  the  E*s, 


3.  As  a provisional  step,  solve  the  global  continuity  equation  for 
the  p’s  and  update  the  pressure  array, 

4.  Solve  the  momentum  equation  for  the  w’s, 

5.  Re-solve  the  global  continuity  equation  for  the  p’s  and  again 
update  the  pressure  array. 


Steps  1-5  constitute  one  iteration;  the  values  of  the  dependent  variables 
at  the  new  time  level  are  taken  to  be  those  obtained  on  the  third  iter- 
ation. Of  the  flow  field  points  investigated,  the  third  iteration  pro- 
duced convergence  to  at  least  five  significant  digits. 


Stability  of  this  integration  could  be  maintained  only  by  observing 
the  well-known  Courant  condition. 

At  < Ax/ (a+u) 

max 

which  is  normally  applicable  to  explicit  techniques-  Furthermore,  the 
cell  Reynolds  number  limitation, 


Re 


cell 


uReAz  < 2 


became  a major  barrier.  All  attempts  to  circumvent  these  two  restrictions 
ended  in  failure.  The  difficulty  stems  primarily  from  the  type  of  pro- 
blem posed  here.  A system  of  coupled  parabolic  partial  differential 
equations,  exhibiting  an  infinite  speed  of  propagation  (all  characteris- 
tic directions  collapse  to  t=constant) , is  solved  between  two  fixed 
boundaries.  The  closed  computational  domain  traps  (and  retains)  numeri- 
cal errors  as  well  as  the  desired  physical  quantitites.  Furthermore, 
rapid  heating  of  the  chamber  excites  nonlinear  acoustic  oscillations  in 
the  gas  which  change  significantly  on  a short  time  scale.  This  combina- 
tion presents  a formidable  challenge  to  any  numerical  method  based  on 
large  time  steps. 
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III.  RESULTS 


h 

Numerical  predictions  for  ignition  behavior  are  presented  for  a 
simple  two-species  system  representing  the  ultimate  in  a premixed  re- 
action. Initially,  the  chamber  contains  only  species  A which  is  capable 
of  the  single,  one-step,  irreversible  reaction, 

A 2B. 

As  would  be  the  case  for  an  elementary  reaction,  the  reaction  rate  is 
assumed  linearly  proportional  to  the  concentration  of  A.  Hence,  the 
^ chemical  production  rate  term  required  in  Eq.  (4)  can  be  written  as 

03^  = -Z*pY^(L*/a^*)  exp[-E*/R^*T*]. 

In  keeping  with  the  simplicity  of  this  example,  the  thermal  conducti- 
vity, viscosity,  diffusion  coefficient,  and  specific  heats  were  all 
held  constant  at  the  representative  values  shown  in  Table  1. 


* 


Table  1.  Parameters  Held  Constant 
30  gm/gm-mole  = 15  gm/gm-mole 

D 

.331  cal/gm-°K  c = .662  cal/gin-°K 

P B 

2.0  X 10  ^ cal/cm-sec-°K 

-1  2 

5.0  X 10  cm  /sec 

5.0  X 10  ^ gm/cm-sec  ; y = 1 
-2 

2.5x10  cm 
3.2238  X 10^  cm/sec 

300  °K 
1965.4 


A parametric  study  was  conducted  to  estimate  the  influence  of^the  three 
important  kinetic  parameters,  Z*  the  pre-exponential  factor,  E the 
activation  energy,  and  AH*  the  heat  of  reaction.  The  values  considered 
are  listed  in  Table  2. 


19 


Table  2.  Results  of  Parametric  Study 
(q^aii  “ cal/cm  -sec) 


* 


Case 

E 

kcal 

mole 

* 

^ -1 
sec 

AH 

(exothermic) 

kcal 

mole 

J^mq 

at  go /no-go 
(n-d) 

Std. 

35 

lO^^ 

35 

127.87 

CVl 

35 

10^^ 

35 

174.7 

CV2 

35 

10^^ 

25 

185.4 

CV3 

35 

10^^ 

30 

179.3 

CV4 

35 

10^3 

30 

131.0 

CVS 

35 

10^" 

30 

97.6 

CV6 

30 

10^^ 

30 

115.8 

CV7 

30 

10^3 

30 

83.0 

The  results  to  follow  are  displayed  in  machine-drawn  plots,  con- 
structed by  connecting  two  points  with  a straight  line  segment.  None 
of  the  curves  were  smoothed  or  fitted  with  polynomials  or  splines.  The 
non-dimensional  axial  distance,  0 ^ z 1,  is  divided  into  50  cells. 

For  reference  purposes.  Fig,  5 shows  the  temperature  distribution 
in  the  chamber  for  the  special  case  when  q^all  held  constant  at 
10  cal/cm^-sec  and  the  reaction  rate  is  set  to  zero.  The  parameter 
values  listed  on  the  right-hand  side  of  the  figure  are  non-dimensional 
times  from  the  onset  of  heating;  one  non-dimensional  time  unit  is  the 
time  required  for  a small  amplitude  disturbance  to  travel  the  length 
of  the  chamber,  assuming  a uniform  temperature  of  300°K  and  = 1. 
Initial  conditions  are  u*  = 0,  p*  = 1 atm.,  and  = 1.  The  results 
show  that  SSO^'K  is  the  maximum  temperature  achieved  during  a heating 
time  of  145  time  units  when  species  A is  inert. 

Figure  6 shows  the  temperature  distribution  in  the  chamber  for  the 
kinetic  parameters  denoted  as  the  standard  case.  The  heating  rate  re- 
mained constant  at  10  cal/cm^-sec.  The  temperature  at  z=0  exceeds  850°K 
at  a time  of  128,  and  the  distribution  shows  evidence  of  the  classic 
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Figure  5.  Temperature  Distribution  for  Standard  Case  (No  Reaction) 
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Figure  6,  Temperature  Distribution  for  Standard  Case 


"bulge”  near  the  heated  wall  as  discussed  by  Merzhanov  and  Averson  • 

This  indicates  the  onset  of  the  reaction.  These  temperature  profiles, 
however,  give  no  indication  of  the  oscillatory  nature  of  the  confined  gas. 
Rapid  heating  of  the  chamber  leads  to  a thermally-driven  convection  pro- 
cess which  creates  the  velocity  oscillation  pattern  shown  in  Fig.  7. 

The  plot  shows  the  development  of  the  wave  over  one-half  of  the  cycle; 
as  indicated  by  the  shape  of  the  profile,  the  frequency  is  close  to 
that  of  the  fundamental.  Note  that  all  velocities  are  positive.  Hence, 
hot  gas  is  continuously  carried  away  from  the  heated  end  of  the  chamber. 
This  is  verified  by  the  development  of  the  density  distribution  shown  in 
Fig.  8.  The  counterpart  of  the  bulge  in  the  temperature  distribution 
near  z=0  is  visible  in  the  density  distribution  at  t=^128.  The  extent  of 
the  chemical  reaction  at  this  time  can  be  visualized  with  the  reactant 
mass  fraction  distribution  displayed  in  Fig.  9.  The  results  indicate 
the  reaction  at  the  boundary  is  approximately  two  percent  complete. 

As  a first  step  toward  establishing  the  critical  go/no-go  ignition 
point,  the  time  t^q  which  begins  the  shut-down  portion  of  the  external 
heat  flux  curve  (see  Fig.  3)  was  set  to  127.0.  The  temperature  distri- 
bution in  Fig.  10  clearly  shows  the  reaction  was  not  self-sustaining, 
and  the  temperature  profiles  decay  between  two  adiabatic  walls.  An  in- 
teresting aspect  of  this  cooling  process  is  the  effect  on  the  velocity 
oscillation  pattern.  The  pressure  near  the  heated  boundary  immediately 
senses  the  falling  temperatures  and  sends  an  expansion  wave  across  the 
chamber.  This  expansion  wave  inverts  the  velocity  oscillation  to  create 
the  pattern  in  Fig.  11,  which  shows  the  development  of  the  distribution 
during  one-half  of  the  oscillation  cycle.  The  negative  velocities  mean 
the  direction  of  the  flow  is  toward  the  heated  boundary;  thus,  cold  gas 
is  brought  in  to  help  "put  out  the  fire."  This  is  again  confirmed  by 
the  density  distribution  displayed  in  Fig.  12.  The  progress  of  the  re- 
action can  be  followed  in  Fig.  13  which  shows  the  reactant  mass  fraction 
distribution  during  the  same  time  interval.  The  reaction  near  the  wall 
continues  to  three  percent  completion  (at  t~133)  even  though  the  tempera- 
tures are  decreasing,  but  then  convection  and  diffusion  begin  to  domi- 
nate. The  result  is  an  increase  in  the  mass  fraction  of  A near  the 
wall,  as  the  product  B is  spread  throughout  the  chamber. 

The  behavior  near  the  critical  ignition  point  can  be  seen  more 
clearly  by  plotting  the  temperature  of  the  gas  at  z=0  as  a function  of 
time,  for  a series  of  runs  varying  only  t^q.  Such  a plot  is  shown  in 
Fig.  14.  The  extreme  sensitivity  of  this  pre-mixed  reaction  to  the 
cut-off  time  of  the  external  stimulus  is  evident;  the  critical  time  of 
127.87  leads  to  a hang-fire  situation  where  the  runaway  reaction  occurs 
36  time  units  after  t^q.  A nearly  identical  phenomenon  was  predicted 

*0 

MerzhanoVjy  A.  G.^  and  Averson^  A,  ^^The  Present  State  of  the 

Thermal  Ignition  Theory:  An  Invited  Review^\  Comb.  Elaine^  Vol.  16^ 
pp  89-124  (1971). 
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Figure  7.  Velocity  Distribution  for  Standard  Case  At  Time  -107 
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Figure  8.  Density  Distribution  for  Standard  Case 
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Figure  9.  Reactant  Mass  Fraction  Distribution  for  Standard  Case 
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Figure  10.  Temperature  Distribution  for  Standard  Case  (t  =127.0) 
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Figure  12.  Density  Distribution  for  Standard  Case  (t  =127.0) 
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Figure  13.  Reactant  Mass  Fraction  Distribution  for  Standard  Case  (t  =127.0) 
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Figure  14.  Temperature-Time  History  At  Heated  Wall  (Varying  t ) for  Standard  Case 


by  Bradley  in  his  i-gnition  study  of  a reactive  solid  with  no  contri- 
bution from  the  gas-phase. 

Evidence  that  the  runaway  reaction  is  self-sustaining  after  the 
long  hang-fire  Ctmq=127.87)  is  given  by  the  temperature  distribution 
plotted  in  Fig.  15  at  a time  of  approximately  164.  The  results  indicate 
a rapid  rise  to  nearly  26G0^K  as  the  flame  forms  at  the  z=0  boundary. 

In  this  short  time  interval,  the  temperature  is  increasing  much  faster 
than  the  density  can  decrease.  Hence,  a significant  change  in  pressure 
should  occur.  This  is  indeed  the  case,  as  shown  in  Fig.  16.  The  rapidly 
increasing  pressure  profiles  coalesce  into  a shock  wave.  Unfortunately, 
the  numerical  integration  method  is  sensitive  to  the  fact  that  continu- 
ous differential  equations  cannot  be  integrated  directly  across  a dis- 
continuity and  this  halts  the  solution.  It  is  concluded,  however,  that 
successful  ignition  of  a pre-mixed  gas  with  the  kinetic  parameters  of 
the  standard  case  will  lead  to  a detonation  wave. 

The  rapid  heat  release  which  drives  the  shock  wave  formation,  dis- 
played in  Figs.  15  and  16,  has  a dramatic  influence  on  the  gas  velocity 
oscillation  pattern  in  the  chamber.  This  is  shown  in  the  multiple 
plots  of  Fig.  17  for  the  standard  case  when  t^q=128.0  (see  Fig.  14  for 
the  wall  temperature- time  history) . Temperature  and  velocity  distri- 
butions are  plotted  in  Figs.  17(a)  and  17(b)  for  the  time  interval  when 
the  runaway  reaction  is  just  beginning.  At  the  time  of  133.14,  the 
velocity  distribution  is  a low  amplitude  oscillation  pattern  similar  to 
that  shown  in  Fig.  7.  However,  the  internal  heat  generation  near  the 
z=0  boundary  creates  an  increasing  temperature  and  pressure  field  which 
halts  the  oscillation  and  forces  a monotonically  swelling  velocity  pro- 
file. This  continues  at  an  accelerated  pace,  as  shown  in  Figs.  17(c) 
and  17 (d) , and  the  maximum  velocity  moves  closer  toward  the  boundary. 

This  will  become  a ”blast-wave"  velocity  profile  which  must  follow  be- 
hind the  forming  detonation  wave. 

A set  of  kinetic  parameters  which  does  not  create  a detonation 
wave  can  be  found  in  Case  CV2.  Compared  to  the  standard  case,  the  pre- 
exponential factor  is  reduced  by  a factor  of  10  and  the  exothermic 
heat  release  is  10  kcal/mole  lower.  This  combination  significantly  in- 
creases the  time  to  ignition  as  shown  in  the  temperature-time  history 
plot  of  Fig.  18.  Once  ignited,  this  case  produces  a much  less  severe 
temperature  distribution  as  shown  in  Fig.  19.  The  flame  temperature  is 
approximately  1800°K,  compared  to  2600®K  for  the  standard  case.  Figure 
20  shows  the  reactant  mass  fraction  distribution  as  the  flame  is  forming. 
An  estimate  of  the  thickness  of  this  pre-mixed  flame  is  about  50  microns. 
This  may  be  reasonable  considering  that  a complex  methane-air  diffusion 
flame  is  600-800  microns.  Figure  21  is  a plot  of  the  pressure  distri- 
bution during  flame  formation.  The  last  profile  is  similar  to  the 


Bradley^  H.  7/.  ^^Theory  of  Ignition  of  a Reactive  Solid  by  Con- 

stant Energy  Flnx*^^  Comb.  Sci.  Tech.:,  Vol.  2^  pip  11-20^  (1970). 
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Figure  15.  Temperature  Distribution  for  Standard  Case  (t  =127.87)  Showing  Flame  Formation 
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Figure  16.  Pressure  Distribution  for  Standard  Case  (t  =127.87)  Showing  Shock  Wave  Formation 
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Figure  17a.  Temperature  Distribution  for  Standard  Case  (t  =128.0)  As  Runaway 
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Figure  17b.  Velocity  Distribution  for  Standard  Case  (t  =128.0)  as  Runaway 
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Figure  17c.  Temperature  Distribution  for  Standard  Case  (t  =128.0)  As  Runaway 
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Figure  17d.  Velocity  Distribution  for  Standard  Case  (t  =128.0)  As  Runaway 
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Figure  18.  Temperature-Time  History  At  Heated  Wall  (Varying  t ) For  Case  CV2 
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Figure  19.  Temperature  Distribution  for  Case  CV2  (t  =185.5)  Showing  Flame  Formation 
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classical  description  of  a laminar  deflagration  wave  in  which  the  pres- 
sure decreases  slightly  across  the  combustion  front.  Of  course,  a com- 
pression wave  must  precede  the  subsonic  deflagration  in  order  to  satis- 
fy the  wall  boundary  condition  at  z=0.  The  slight  ”raggedness’*  of  the 
predicted  pressure  in  this  region  is  non-physical,  and  is  the  direct 
result  of  violating  the  cell  Reynolds  number  limitation  discussed 
earlier.  Increased  grid  point  resolution  would  eliminate  the  problem. 

Correlation  of  Results 


The  results  obtained  in  the  parametric  study  for  time  to  ignition 
are  summarized  in  Table  2.  The  critical  value  of  ignition  time  for  each 
case  is  taken  to  be  the  value  of  t^^q  which  yields  the  longest  hang-fire. 
The  temperature-time  history  for  each  case  can  be  found  in  Appendix  A. 

A better  perspective  of  the  sensitivity  of  ignition  time  to  a change  in 
one  of  the  kinetic  parameters  can  be  seen  with  a simple  power  law  corre- 
lation. These  correlations  are  given  below: 


(a) 


(b) 


(c) 


ignition  time  vs. 
t 


* 2 9 

a (E 


mq 

ignition  time  vs. 

t a (AH  ) 
mq 

ignition  time  vs. 

t a (Z  ) 
mq 


activation  energy 


heat  of  reaction 


pre-exponential  factor 


Thus,  an  increase  in  activation  energy  will  substantially  increase  the 
time  to  ignition,  while  an  increase  in  the  heat  of  reaction  and/or  the 
pre-exponential  factor  will  decrease  the  time  to  ignition. 


In  ignition  studies  of  all  types,  the  question  most  frequently  ad- 
dressed concerns  the  dependence  of  ignition  time  on  external  heat  flux. 

To  answer  this  question  for  a confined  pre-mixed  gas,  the  present  investi- 
gation conducted  a special  series  of  runs  based  on  the  standard  case. 

(The  results  for  temperature-time  history  are  included  in  Appendix  B) . 

The  constant  portion  of  q^rall^^^  varied  in  small  increments  over  the 
limited  range  of  10-25  cal/cm^-sec.  The  predictions  for  t at  go/no- 
go  as  a function  of  the  constant  portion  of  correlated  exactly  as 

the  power  law. 


t 

mq 


wall 


-m 


m=1.89 
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An  interesting  comparison  can  be  made  with  the  results  of  two  different 
solid  propellant  ignition  theories.  Bradley  and  Williams^^  postulated  a 
heterogeneous  reaction  between  a gas-phase  oxidizer  and  the  surface  of 
the  solid  propellant.  Their  correlation  produces  the  above  result  with 
m=1.86,  for  the  activation  energy  used  in  the  standard  case  of  this 
paper.  Kashiwagi^  postulated  a two  reaction  model;  the  solid  fuel 
pyrolyzes  (endothermically)  at  the  surface  and  then  reacts  in  a gas- 
phase  diffusion  flame  with  the  ambient  oxidizer.  The  result  of  his 
correlation  produces  the  exponent  m=1.9,  in  the  heat  flux  range  quoted 
above. 

All  three  correlations  yield  nearly  the  same  value  for  the  exponent. 
In  fact,  this  value  is  close  to  2.00  which  is  the  value  obtained  from  the 
well-known  solution  to  the  heat  conduction  equation  for  an  inert  solid 
(e.g.,  see  discussion  in  Ref.  3).  The  important  implication  here  is 
that  a plot  of  ignition  time  versus  external  heat  flux  cannot  be  used 
to  distinguish  the  mechanism  controlling  ignition. 

IV.  CONCLUSIONS 

1.  Numerical  predictions  of  the  gas-phase  ignition  sequence  seem  , 
qualitatively  correct.  In  general,  the  results  show  that: 

a.  A confined  premixed  gas  capable  of  a single,  one-step, 
irreversible,  chemical  reaction  can  sustain  a hang-fire  ignition,  i.e., 
the  runaway  reaction  occurs  many  time  units  after  the  external  heat 
source  (ignition  stimulus)  is  terminated. 

b.  A hang-fire  ignition  forms  the  boundary  separating  *'go"  and 


c.  The  acoustic  response  of  the  chamber  plays  a definite  role 
in  the  ignition  of  a confined  reactive  gas. 

d.  Successful  ignition  of  the  confined  premixed  gas  leads  to 

a thin  flame  zone,  assuming  that  all  transport  processes  are  laminar 
and  that  the  controlling  chemical  reaction  is  a single,  one-step, 
irreversible  decomposition.  When  the  exothermic  heat  of  reaction  is 
sufficiently  large,  the  numerical  solution  predicts  the  formation  of 
a detonation  wave. 


Bvadtey^  H.  and  VlilZioms^  F.  A,  ^ '^Theory  of  Radiant  and  Ey- 

pevgolic  Ignition  of  Solid  Propellants^^ ^ Comb.  Soi.  Tech.:,  Vol.  2^, 
pp  41-52  (1970). 
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2.  Correlations  of  the  numerical  results  were  established  to  show 
that: 

* 

a.  An  increase  in  activation  energy  (E  ) will  substantially  in- 
crease the  time  to  ignition. 

* 

b.  An  increase  in  the  exothermic  heat  of  reaction  (AH  ) will  de- 
crease the  time  to  ignition. 

* 

c.  An  increase  in  the  pre-exponential  factor  (Z  ) will  decrease 
the  time  to  ignition. 

3.  The  dependence  of  ignition  time  on  external  heat  flux  cannot  be 
used  to  distinguish  a simple  premixed  gas-phase  reaction  from  other 
mechanisms  which  may  control  ignition. 
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APPENDIX  A 


Ignition  Results  for  Parametric  Study 


V' 


This  appendix  contains  the  numerical  results  obtained  near  the 
critical  ignition  point  for  each  case  considered  in  the  parametric  study 
(for  the  Standard  Case,  see  Fig.  14,  page  31).  Each  of  the  figures  to 
follow  pertains  to  a single  case  and  displays  multiple  curves  of  gas 
temperature  adjacent  to  the  heated  wall  (2*0)  as  a function  of  time; 
the  only  parameter  varied  is  the  cut-off  time,  of  the  external 

heat  flux.  The  constant  portion  of  the  external  heat  flux,  q^all^^^* 
is  held  constant  at  10  cal/cm^-sec.  The  value  of  which  yields  the 

and  was  recorded 


longest  hang-fire  is  taken  as  the  ”time  to  ignition' 
in  Table  2 on  page  20. 


mq 
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Figure  A-1  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  1 
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Figure  A-2  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  2 
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Figure  A-3  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  3 
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Figure  A~4  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  4 
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Figure  A-5  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  5 
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Figure  A-6  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  6 
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Figure  A-7  Temperature-Time  History  at  Heated  Wall  (Varying 

for  Case  CV  7 


APPENDIX  B 


Ignition  Results  for  Variable  External  Heat  Flux 


This  appendix  contains  the  numerical  results  obtained  near  the 
critical  ignition  point  for  the  study  of  ignition  time  dependence  on 
external  heat  flux  level.  All  of  the  figures  to  follow  are  based  on 
the  set  of  parameters  denoted  as  the  Standard  Case.  Each  figure 
corresponds  to  a particular  value  for  the  constant  portion  of  the  ex- 
ternal heat  flux,  q^all^^^  (see  Fig.  3,  page  16).  The  value  of  this 
constant  portion  is  varied  over  the  range  10-25  cal/cm^-sec.  Similar 
to  Appendix  A,  the  results  are  displayed  as  multiple  curves  of  gas 
temperature  adjacent  to  the  heated  wall  (z=0)  as  a function  of  time; 
the  only  parameter  varied  is  the  cut-off  time,  t^^^,  of  the  external  heat 
flux.  The  value  of  which  yields  the  longest  hang-fire  is  taken 
as  the  "time  to  ignition"  and  was  used  to  obtain  the  correlation, 
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Figure  B-1  Temperature-Time  History  at  Headed  Wall  (Varying 
for  "^12.5  cal/cm  -sec 
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Figure  B-2  Temperature-Time  History  at  Headed  Wall  (Varying 

for  ~ 15.0  cal/cm  -sec 
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Figure  B-3  Temperature-Time  History  at  Headed  Wall  (Varying 
for  ~ 20.0  cal/cm  -sec 


1 100 


|-LiJI[LUJ[K<|-DKUJ  ^ 


61 


Figure  B-4  Temperature-Time  History  at  Heated  Wall  (Varying 

for  “ 25.0  cal/cm  -sec 
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